Real time evolution at finite temperatures with operator space matrix product states 
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We propose a method to simulate the real time evolution of one dimensional quantum many- 
body systems at finite temperature by expressing both the density matrices and the observables as 
matrix product states. This allows the calculation of expectation values and correlation functions 
as scalar products in operator space. The simulations of density matrices in inverse temperature 
and the local operators in the Heisenberg picture are independent and result in a grid of expectation 
values for all intermediate temperatures and times. Simulations can be performed using purely real 
arithmetics with only polynomial growth of computational resources in inverse temperature and 
time for integrable systems. The method is illustrated for the XXZ model and the single impurity 
Anderson model. 

PACS numbers: 05.10.-a, 07.05.Tp, 03.67.Bg, 78.47.-p, 75.10. Jm, 75.10.Pq 



The simulation of strongly interacting quantum many- 
body systems is, in general, a computationally difficult 
task. Even in the case of one-dimensional quantum sys- 
tems with finite range interactions, where ground state 
properties have become readily accessible using the den- 
sity matrix renormalization group (DMRG) algorithm 
PP, accessing the dynamics remains hard. While, the 
matrix product states (MPS) [2r[4] that are at the core 
of DMRG manage to describe the ground states due to 
their weak entanglement [5], the entanglement increases 
in the real time evolution [51 [SHI] making the study of dy- 
namics computationally challenging even for pure states. 

Even harder is dynamics at finite temperatures. MPS 
based methods operate on the level of wave functions. 
Systems at thermal equilibrium, described by a mixture 
of quantum states, thus need to be purified by being de- 
scribed as a pure state in an enlarged system, consisting 
of the system and an auxiliary environment. The simula- 
tion is then performed on this purified system (4[ [8] and 
by tracing over the degrees of freedom of the environ- 
ment, the properties of the original system are obtained. 
This idea was first implemented in Ref. [H |H| for den- 
sity matrices of thermal quantum systems and systems 
far from equilibrium. The same concept was used in the 
simulation of dynamical properties of thermal systems 
where the real time evolution was computed on the level 
of the purified system and the density operator for the 
original system was obtained by tracing over the environ- 
ment main]. 

A disadvantage of these purification based approaches 
is that they require additional information on the en- 
vironment, and especially on its dynamics. In gen- 
eral, these methods are accompanied by an exponential 
growth of computational resources with time and elabo- 
rate procedures are required to control the fast growth 
of entanglement. Recent advances in the application to 
the XXZ model have shown that a smart choice of the 
environmental dynamics can significantly extend the ac- 



cessible simulation times [TTJ Q2] . However, the question 
of the optimal environment and its dynamics is still open. 
Even finding a good choice of the environmental dynam- 
ics requires deep knowledge of the system and is far from 
being straightforward. 

Here we propose an alternative strategy, based on oper- 
ator space concepts and matrix product states [13] . Simi- 
larly to the purification approach, we describe the density 
matrices as many-body states; however, our states are 
not superpositions of quantum states but superpositions 
of operators [HI 03] ■ In the same way we also describe the 
observables and their dynamics, known as the time de- 
pendent DMRG in the Heisenberg picture 13, 15 . The 
key feature of our approach is the combination of both 
concepts where the expectation value of an observable at 
a given time and temperature is calculated as a scalar 
product of two vectors in the operator space, one repre- 
senting the density matrix and the other the real time 
evolution of an operator. This way, the simulation of 
density matrices and the real time evolution of opera- 
tors are decoupled and can be performed independently 
such that only two simulations are sufficient to calculate 
a grid of expectation values at various intermediate (in- 
verse) temperatures and times. 

The computational costs to simulate the density ma- 
trices grow polynomially in time due to a low degree of 
correlation at thermal equilibrium |14j , although such de- 
scription has mostly been used to simulate systems far 
from the equilibrium [16] . Similarly, the simulation of lo- 
cal operators also exhibits polynomial growth of resources 
[13j for integrable systems, corresponding to logarithmic 
growth of the operator space equivalent of the entan- 
glement entropy (OSEE) [T7], unlike the simulation of 
quantum pure states where the entanglement entropy in 
the time evolution typically grows linearly [TH] (except 
for e.g. local quenches). Our method can therefore be 
used to simulate the real time evolution of local operators 
in integrable quantum systems in polynomial time. 
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Method - Let us consider a quantum system of n sites 
with a Hilbert space T-L. The expectation value of an op- 
erator a : H — > H in a mixed state described by a density 
matrix p, is given by (a) = tr(pa) with tr(p) = 1. Like 
the operator a, the density matrix is also a linear map 
over T-L and these maps form the operator space algebra 
JC with an inner product (a\b) = (dimH)-Hx(a H b). The 
expectation value of a can thus be given as 

(a) = (p\a)/(p\e), (1) 

where the bracket notation \a) is used when referring to 
a G JC and |e) corresponds to the fully mixed state (the 
identity). In thermal equilibrium for inverse tempera- 
ture /? = (fceT) -1 , the expectation values are given by 
using \p{(i)) = |e _/3ff ) as the density matrix whereas the 
time evolution of Eq. ([!]) is obtained by replacing a by 
the Heisenberg operator a(t) = e ltH ae~ ltH . Clearly, the 
simulations for p(/3) and a(t) can be done independently. 

The density matrices p{j3) are obtained by the imagi- 
nary time evolution, starting from the infinite tempera- 
ture equilibrium state |e). By defining a map over the 
operator space, x : /C — > /C, which left-multiplies an 
operator by H, i.e. x\ a ) = \Ha), the density matrix 
|p(/3)) = \e~P H ) becomes 

|p(0)) = e-«|e). (2) 

The thermal super-operator \ nas the same form and the 
same range of interactions as H , if the basis for JC is local. 

Similarly, the Heisenberg evolution a(t) = e ltH ae~ ltH 
can be simulated by introducing a linear map defined 
as H\x) = \[x,H]) for x € JC. It was shown [17] that 
the Heisenberg evolution is hence transformed to the 
Schrodinger evolution in the operator space 

\a(t)) = e- it6 \a(0)). (3) 

Like x, the super-operator H has the same form as H. 

Let us now define a suitable basis for JC. While any ba- 
sis could in principle be used, it is advantageous to choose 
a basis given by direct products of local basis operators 
{Pvj } pertaining to the site j , as P„ = p\}} <gipl 2 J <g> • • • ®p^} . 
These operators are chosen such that they form an or- 
thonormal set. A suitable choice for {pv}} depends 
on the type of simulation and a different choice can 
be made for the thermal |2]) and the real time evolu- 
tion ([3]). For the former, the map x, represented by 
[x]ju,w = (dim - H)~ 1 tr(P^x-Pf)7 is real if {P u } are real 
which in turn allows us to simulate the thermal density 
matrices ([2| with real arithmetics. 

For the real time evolution (|3| , where H is represented 
by a matrix [H]^ = (dim'H)"Hr([P / f , P V ]H), a different 
choice for {P„} is made to maintain the formulation in 
real arithmetics. If P^ — P M , then [H\n tV = —[H] v>li and 
H is represented by a hermitian skew-symmetric matrix. 
Together with the imaginary unit in the exponent of (p3J) , 



the Heisenberg evolution is generated by a real valued 
map. Finally, due to different operator space bases used 
in the thermal and the Heisenberg evolution, the calcu- 
lation of expectation values must be preceded by a basis 
transformation, which in our case is a product operator 
and has no significant effects on the computational costs. 

For reference we list a few typical operator space bases. 
In the case of spin-1/2 models, spinless fermions or hard 
core bosons the basis is given by the products of Pauli 
matrices pjfj = a^ 3 for pj 6 {0,1,2,3}. This basis is 
hermitian and thus suitable for the real time evolution 
whereas the thermal evolution is done using a real lo- 
cal basis {a , a 1 , —ia 2 , a 3 } . The local basis for spin-1 
systems is provided by Gell-Mann matrices whereas for 
spin-3/2 it can be chosen as {S^ 11 ^ = a a ® a b }. In all 
cases, the local basis is hermitian and its real represen- 
tation can be found in a straightforward way. 

Let us now focus on the MPS ansatz for the operators. 
For simplicity we shall assume a spin-1/2 model although 
the results are readily generalized to any many-body 
quantum system with a local dimension d. The oper- 
ators P Vl ,..., v „ — cr" 1 ■ ■ ■ a^™ form the basis set \v) = \P V ) 
and a € JC is given in terms of a MPS [T3] 

\a)= trtAW^.-AW^K...,^), 

L>1 , . .. ,1/ n 

with real matrices of dimension at most D x D 

(see e.g. [1]) where D is the bond dimension. The evolu- 
tions © and (|]) are simulated using the Suzuki- Trotter 
decomposition and local updates of the MPS, known 
as the time dependent DMRG algorithm [5HZ] which is 
particularly efficient when H contains at most nearest- 
neighbor interactions. 

The thermal expectation value (A(t))^ given by Q is 
obtained by combining the independent simulations ([2| 
and |3]). Using the same principles, we can simulate 
(a{t)b)p where a(t) and b are MPS, by constructing a 
matrix product operator B : x 1— > bx. The result is given 
by (ba(t))p = {p(0)\B\a(t))/(p(f3)\e) which in the MPS 
language corresponds to calculating an expectation value. 
Similarly (a(t)b)p is obtained. 

Results - The method proposed here can be used 
to simulate the real time evolution of a wide range of 
models at finite temperature, for local operators such 
as the spin projection at a given site, a local current, 
or correlations of local operators (e.g. Green's func- 
tions). We start by analyzing the complexity of both 
constituents of the method, the evolution of the den- 
sity matrices and the time evolution of operators in the 
Heisenberg picture. We quantify the complexity in terms 
of the OSEE [17] which is an operator space analogue of 
the entanglement entropy for quantum states, defined as 
S*(x) = -tr[i?log 2 7?] for R = tr B (|jc)(a:|). For density 
matrices, it was observed [2] that the OSEE is related 
to the quantum mutual information which quantifies the 
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FIG. 1. Operator space entanglement entropy (OSEE) in the 
evolution of density matrices p(/3) in the inverse temperature 
ft (top) and the local current operator j n /i{t) (bottom) for 
the XXZ model on 100 sites with A indicated in the legends. 
Bond dimensions are D = 1024 and D = 1200, respectively. 
A logarithmic dependence in (imaginary) time is observed in 
both cases (see the insets in a semi- log scale). 

total correlations between two parts of the system. 

We first consider a spin- 1/2 XXZ model with open 
boundary conditions described by the Hamiltonian 

n-i 

H = E (*f °f+i + °? ff f+i + Aofof+i) > ( 4 ) 
i=i 

where A is the anisotropy parameter. The OSEE of the 
initial state p(j3 = 0) is zero (the fully mixed state is sep- 
arable). As shown in Fig. [T] (top), we observe logarith- 
mic growth of the OSEE in the inverse temperature [TJ] . 
Furthermore, we simulate the Heisenberg evolution ([3| 
of the spin current flowing between sites m and 77?, + 1 

jm = «Oro°m+l - °"m CT m+l)> with a m = (°"m ± ^m)/ 2 

for m = n/2. The results in Fig. [T] (bottom) show the 
logarithmic growth of OSEE which is an extension of the 
results in [13 to non-linear integrable systems. Simi- 
lar to the evolution of density matrices, we observe an 
increase of OSEE for larger anisotropics A where the 
simulation is more expensive. The logarithmic growth 
of OSEE corresponds to polynomial growth of the bond 
dimension |131 117) and the simulation can be performed 



FIG. 2. Current correlation function for the XXZ model on 
200 sites for various A and fixed /3 = 0.25. Different symbols 
correspond to maximal bond dimensions D = 800, 1000, 1200 
(real time). The error in p{0) is negligible for D = 200. 

efficiently. We note that the computational cost is ex- 
pected to grow exponentially |13| for non-integrable sys- 
tems, like in other similar methods. 

The simulation of the density matrices and the local 
operators can now be used to calculate e.g. the spin 
current correlation function (jj n /2('t)) [}■ The results are 
shown in Fig. [2j an inverse temperature /3 = 0.25. The 
infinite-time limit of these current-correlations yields the 
Drude weight which can be used to signal ballistic trans- 
port in various integrable models, see [19]. For the XXZ 
chain similar results were obtained in Refs. Q21 [501 l2"Tj . 
a detailed discussion is however beyond the scope of this 
manuscript and will be presented elsewhere. When com- 
pared to Ref. [12] (with our time scale multiplied by a 
factor of 4 due to the spin operators) we observe that 
longer times can be reached. This appears particularly 
pronounced for large values of A with reduced maximal 
simulation times as shown in Fig. [2] 

As another illustration of the method we calculate the 
Green's functions of the single impurity Anderson model 
which can be mapped [22] [23] to an XXZ-like model 

oo 

H = ^ r j(°/°j+i + + ^"-o^i + e/(no + n{) 

j=-oo 

(5) 

with rij = a- (Tj and the interaction strength U. The 
hopping parameters r_i = n are given by the hybridiza- 
tion whereas the rest of Tj are determined from the dis- 
cretization of the bath (and To = 0), see [23] for de- 
tails. The interaction term is only present between sites 
and 1 (corresponding to spin-f and spin- 4- of the im- 
purity) which allows longer acessible times and lower 
temperatures. We simulate the single particle Green's 
function G(t) = —i({fl,ff(t)})p where is the anni- 
hilation operator at the impurity for spin-f, written as 
ft = (n.j<a a j) <J o d ue to Jordan- Wigner transforma- 
tion. Despite the nonlocality, the Heisenberg evolution 
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and SIMCOFE and by SFB ViCoM (FWF project ID 
F4104). Simulations were run on the Brutus cluster at 
ETH Zurich and on the Vienna Scientific Cluster. 



FIG. 3. Temporal Green's function for the single impurity 
Anderson model |5| for U = 1, Tj = |, To = 0, and the 
particle-hole symmetric e/ = —U/2 where local Coulomb cor- 
relations are most pronounced. 



of this operator can be simulated efficiently [13j \T7\ . We 
actually simulate Majorana operators w = /-]- + /j and 

w' = i(/f — /|) which are hermitian and allow the simula- 
tion in the real arithmetics. The results for the imaginary 
part of the Green's function, shown in Fig. [3] can be used 
to calculate the spectral functions (see e.g. Ref. [2~4"H2l)] ) 
by means of a Fourier transform which will be in detail 
presented elsewhere. The method can also be used to 
study nonequilibrium impurity physics [27], in particu- 
lar current-voltage characteristics of strongly correlated 
nanostructures in the non-linear response regime. 

Conclusions - We have proposed a tensor network 
method to simulate the real time dynamics of one di- 
mensional quantum many-body systems at finite tem- 
peratures. The main advantage of the method is that 
it decouples the real time evolution from the imaginary 
evolution in the inverse temperature such that two in- 
dependent simulations give the results for all combina- 
tions temperature-time. Furthermore, the simulations 
are done in real arithmetics. The method is most useful 
for integrable systems where the costs grow polynomially, 
and systems close to the integrability where the asymp- 
totic exponential growth only appears at later times. The 
proposed framework can be extended in a straight for- 
ward way to two dimensional systems, described in terms 
of projected entangled pair states [2"5] . 
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